In [3]:
# STAT 415/615 Regression (M. Baron)

# Python Lab 7: Matrix Operations and Multivariate Linear Regression

# Define a matrix by entering its elements manually.

import numpy as np

B = np.array([
[1, 4, 7],
[2, 5, 8],
[3, 6, 9]
])

B
Out[3]:
array([[1, 4, 7],
       [2, 5, 8],
       [3, 6, 9]])
In [4]:
# B**2 is an element-wise operation

B**2
Out[4]:
array([[ 1, 16, 49],
       [ 4, 25, 64],
       [ 9, 36, 81]])
In [5]:
# Matrix multiplication, B^2 = B * B

B @ B
Out[5]:
array([[ 30,  66, 102],
       [ 36,  81, 126],
       [ 42,  96, 150]])
In [6]:
# Transposed matrix

B.T
Out[6]:
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
In [7]:
# Joining two matrices side by side (as columns)

np.column_stack((B, B))
Out[7]:
array([[1, 4, 7, 1, 4, 7],
       [2, 5, 8, 2, 5, 8],
       [3, 6, 9, 3, 6, 9]])
In [8]:
# Joining two matrices below each other (as rows)

np.row_stack((B, B))
Out[8]:
array([[1, 4, 7],
       [2, 5, 8],
       [3, 6, 9],
       [1, 4, 7],
       [2, 5, 8],
       [3, 6, 9]])
In [9]:
# Sub-matrix, a part of matrix B

B[0:2, 0:3]
Out[9]:
array([[1, 4, 7],
       [2, 5, 8]])
In [10]:
# Inverting matrices is available in NumPy

np.linalg.inv(B)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[10], line 3
      1 # Inverting matrices is available in NumPy
----> 3 np.linalg.inv(B)

File ~\AppData\Local\anaconda3\Lib\site-packages\numpy\linalg\linalg.py:561, in inv(a)
    559 signature = 'D->D' if isComplexType(t) else 'd->d'
    560 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 561 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    562 return wrap(ainv.astype(result_t, copy=False))

File ~\AppData\Local\anaconda3\Lib\site-packages\numpy\linalg\linalg.py:112, in _raise_linalgerror_singular(err, flag)
    111 def _raise_linalgerror_singular(err, flag):
--> 112     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix
In [11]:
# What happened? Singular matrix means that there is a linear dependence among columns (and among rows)
# of matrix B. Such matrices are not invertible, and they have a determinant equal to det(B) = 0.
# So, we did not make an error, except that we tried to invert a non-invertible matrix.
# Okay, let's change one element, to make sure the inverse matrix exists.

np.linalg.det(B)

B[0, 0] = 100

B
Out[11]:
array([[100,   4,   7],
       [  2,   5,   8],
       [  3,   6,   9]])
In [12]:
# We changed the matrix by adding a "ridge", and now the inverse B^(-1) exists.

np.linalg.inv(B)
Out[12]:
array([[ 0.01010101, -0.02020202,  0.01010101],
       [-0.02020202, -2.95959596,  2.64646465],
       [ 0.01010101,  1.97979798, -1.65656566]])
In [14]:
# Define a matrix from the "mtcars" data set and build a regression model
# that predicts miles per gallon based on the number of cylinders,
# horsepower, axel ratio, weight, and acceleration time.

import statsmodels.api as sm

mtcars = sm.datasets.get_rdataset("mtcars").data

mtcars.head()
Out[14]:
mpg cyl disp hp drat wt qsec vs am gear carb
rownames
Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
In [15]:
# X-matrix of predictors

X = mtcars[["cyl", "hp", "drat", "wt", "qsec"]].values
In [20]:
X[:6]
Out[20]:
array([[  6.   , 110.   ,   3.9  ,   2.62 ,  16.46 ],
       [  6.   , 110.   ,   3.9  ,   2.875,  17.02 ],
       [  4.   ,  93.   ,   3.85 ,   2.32 ,  18.61 ],
       [  6.   , 110.   ,   3.08 ,   3.215,  19.44 ],
       [  8.   , 175.   ,   3.15 ,   3.44 ,  17.02 ],
       [  6.   , 105.   ,   2.76 ,   3.46 ,  20.22 ]])
In [21]:
# Vector of responses

Y = mtcars["mpg"].values

n = len(Y)

# We also need a vector of 1s to include the intercept

one = np.ones((n, 1))

X = np.column_stack((one, X))

# This is matrix X'X

X.T @ X
Out[21]:
array([[3.20000000e+01, 1.98000000e+02, 4.69400000e+03, 1.15090000e+02,
        1.02952000e+02, 5.71160000e+02],
       [1.98000000e+02, 1.32400000e+03, 3.22040000e+04, 6.91400000e+02,
        6.79404000e+02, 3.47556000e+03],
       [4.69400000e+03, 3.22040000e+04, 8.34278000e+05, 1.63722800e+04,
        1.64717440e+04, 8.10921600e+04],
       [1.15090000e+02, 6.91400000e+02, 1.63722800e+04, 4.22790700e+02,
        3.58718960e+02, 2.05691400e+03],
       [1.02952000e+02, 6.79404000e+02, 1.64717440e+04, 3.58718960e+02,
        3.60901070e+02, 1.82809458e+03],
       [5.71160000e+02, 3.47556000e+03, 8.10921600e+04, 2.05691400e+03,
        1.82809458e+03, 1.02934802e+04]])
In [22]:
# Slope β = (X'X)^(-1) X'Y

slope = np.linalg.inv(X.T @ X) @ X.T @ Y

slope
Out[22]:
array([ 2.59452057e+01, -4.89672771e-01, -1.53892029e-02,  1.13077207e+00,
       -3.38279350e+00,  3.50112654e-01])
In [23]:
# We can certainly get the same slopes by the usual regression command "OLS"

X_reg = sm.add_constant(
mtcars[["cyl", "hp", "drat", "wt", "qsec"]]
)

reg = sm.OLS(mtcars["mpg"], X_reg).fit()

print(reg.params)

const    25.945206
cyl      -0.489673
hp       -0.015389
drat      1.130772
wt       -3.382794
qsec      0.350113
dtype: float64
In [24]:
# Our estimated regression equation is

# mpg = 25.95 - 0.49 cyl - 0.015 hp + 1.13 drat - 3.38 wt + 0.35 qsec + ε